home *** CD-ROM | disk | FTP | other *** search
/ Super PC 34 / Super PC 34 (Shareware).iso / spc / UTIL / DJGPP2 / V2 / DJTST200.ZIP / tests / libc / ansi / math / elefunt / tsin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-01  |  7.9 KB  |  332 lines
  1. /* -*-C-*- tsin.c */
  2.  
  3. #include "elefunt.h"
  4.  
  5. /*#     PROGRAM to test sin/cos
  6. #
  7. #     data required
  8. #
  9. #        none
  10. #
  11. #     subprograms required from this package
  12. #
  13. #        machar - an environmental inquiry program providing
  14. #                 information on the floating-point arithmetic
  15. #                 system.  note that the call to machar can
  16. #                 be deleted provided the following five
  17. #                 parameters are assigned the values indicated
  18. #
  19. #                 ibeta  - the radix of the floating-point system
  20. #                 it     - the number of base-ibeta digits in the
  21. #                          significand of a floating-point number
  22. #                 minexp - the largest in magnitude negative
  23. #                          integer such that  float(ibeta)**minexp
  24. #                          is a positive floating-point number
  25. #                 eps    - the smallest positive floating-point
  26. #                          number such that 1.0+eps != 1.0
  27. #                 epsneg - the smallest positive floating-point
  28. #                          number such that 1.0-epsneg != 1.0
  29. #
  30. #        ran(k) - a function subprogram returning random real
  31. #                 numbers uniformly distributed over (0,1)
  32. #
  33. #
  34. #     standard fortran subprograms required
  35. #
  36. #         abs, alog, amax1, cos, float, sin, sqrt
  37. #
  38. #
  39. #     latest revision - december 6, 1979
  40. #
  41. #     author - w. j. cody
  42. #              argonne national laboratory
  43. #
  44. #*/
  45. #if 0
  46. #ifdef LDOUBLE
  47. static long double _sinl(long double x)
  48. {
  49.     static long double PIO4L   = 7.8539816339744830961566E-1L;
  50.     static long double DP1 = 7.853981554508209228515625E-1L;
  51.     static long double DP2 = 7.946627356147928367136046290398E-9L;
  52.     static long double DP3 = 3.061616997868382943065164830688E-17L;
  53.     static long double lossth = 5.49755813888e11L; /* 2^39 */
  54.     long double y, z;
  55.     int sign = 1;
  56.     long j;
  57.  
  58.     if( x < 0 )
  59.     {
  60.     x = -x;
  61.     sign = -1;
  62.     }
  63.     if( x > lossth )
  64.     {
  65.     return _MATH_NANL;
  66.     }
  67.  
  68.     y = floorl( x/PIO4L ); /* integer part of x/PIO4 */
  69.     z = y;
  70.     j = z; /* convert to integer for tests on the phase angle */
  71.     /* map zeros to origin */
  72.     if( j & 1 )
  73.     {
  74.     j += 1;
  75.     y += 1.0;
  76.     }
  77.     j = j & 07; /* octant modulo 360 degrees */
  78.     /* reflect in x axis */
  79.     if( j > 3)
  80.     {
  81.     sign = -sign;
  82.     j -= 4;
  83.     }
  84.     /* Extended precision modular arithmetic */
  85.     z = ((x - y * DP1) - y * DP2) - y * DP3;
  86.  
  87.     if( (j==1) || (j==2) )
  88.     {
  89.     y = cosl(z);
  90.     }
  91.     else
  92.     {
  93.     y = sinl(z);
  94.     }
  95.  
  96.     if(sign < 0)
  97.     y = -y;
  98.     return(y);
  99. }
  100. #undef sin
  101. #define sin _sinl
  102. #endif
  103. #endif
  104.  
  105. void
  106. tsin()
  107. {
  108.     int i,
  109.     ibeta,
  110.         iexp,
  111.         irnd,
  112.         it,
  113.         i1,
  114.         j,
  115.         k1,
  116.         k2,
  117.         k3,
  118.         machep,
  119.         maxexp,
  120.         minexp,
  121.         n,
  122.         negep,
  123.         ngrd;
  124.  
  125.     float a,
  126.         ait,
  127.         albeta,
  128.         b,
  129.         beta,
  130.         betap,
  131.         c,
  132.         del,
  133.         eps,
  134.         epsneg,
  135.         expon,
  136.         r6,
  137.         r7,
  138.         three,
  139.         w,
  140.         x,
  141.         xl,
  142.         xmax,
  143.         xmin,
  144.         xn,
  145.         x1,
  146.         y,
  147.         z,
  148.         zz;
  149.     machar(&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp,
  150.         &maxexp, &eps, &epsneg, &xmin, &xmax);
  151.     beta = (float) ibeta;
  152.     albeta = ALOG(beta);
  153.     ait = (float) it;
  154.     three = 3.0e0L;
  155.     a = ZERO;
  156.     b = 1.570796327e0L;
  157.     c = b;
  158.     n = 2000;
  159.     xn = (float) n;
  160.     i1 = 0;
  161.  
  162.     /* random argument accuracy tests */
  163.  
  164.     for (j = 1; j <= 3; ++j)
  165.     {
  166.     k1 = 0;
  167.     k3 = 0;
  168.     x1 = ZERO;
  169.     r6 = ZERO;
  170.     r7 = ZERO;
  171.     del = (b - a) / xn;
  172.     xl = a;
  173.  
  174.     for (i = 1; i <= n; ++i)
  175.     {
  176.         x = del * ran(i1) + xl;
  177.         y = x / three;
  178.         y += x;
  179.         y -= x;
  180.         x = three * y;
  181.         if (j == 3)
  182.         {
  183.         z = cos(x);
  184.         zz = cos(y);
  185.         w = ONE;
  186.         if (z != ZERO)
  187.             w = (z + zz * (three - 4.0e0L * zz * zz)) / z;
  188.         }
  189.         else
  190.         {
  191.         z = sin(x);
  192.         zz = sin(y);
  193.         w = ONE;
  194.         if (z != ZERO)
  195.             w = (z - zz * (three - 4.0e0L * zz * zz)) / z;
  196.         }
  197.         if (w > ZERO)
  198.         k1 = k1 + 1;
  199.         if (w < ZERO)
  200.         k3 = k3 + 1;
  201.         w = ABS(w);
  202.         if (w > r6)
  203.         {
  204.         r6 = w;
  205.         x1 = x;
  206.         }
  207.         r7 = r7 + w * w;
  208.         xl = xl + del;
  209.     }
  210.  
  211.     k2 = n - k3 - k1;
  212.     r7 = sqrt(r7 / xn);
  213.     if (j == 3)
  214.     {
  215.         printf("\fTEST OF COS(X) VS 4*COS(X/3)**3-3*COS(X/3)\n\n");
  216.         printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n",n);
  217.         printf("      (" F15P4E "," F15P4E ")\n\n\n", a, b);
  218.         printf(" COS(X) WAS LARGER%6d TIMES,\n", k1);
  219.         printf("            AGREED%6d TIMES, AND\n", k2);
  220.         printf("       WAS SMALLER%6d TIMES.\n\n", k3);
  221.     }
  222.     else
  223.     {
  224.         printf("\fTEST OF SIN(X) VS 3*SIN(X/3)-4*SIN(X/3)**3\n\n");
  225.         printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n", n);
  226.         printf("      (" F15P4E "," F15P4E ")\n\n\n", a, b);
  227.         printf(" SIN(X) WAS LARGER%6d TIMES,\n", k1);
  228.         printf("            AGREED%6d TIMES, AND\n", k2);
  229.         printf("       WAS SMALLER%6d TIMES.\n\n", k3);
  230.     }
  231.     printf(
  232. " THERE ARE%4d BASE%4d SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER\n\n",
  233.         it, ibeta);
  234.     w = -999.0e0;
  235.     if (r6 != ZERO)
  236.         w = ALOG(ABS(r6)) / albeta;
  237.     printf(" THE MAXIMUM RELATIVE ERROR OF" F15P4E " = %4d **" F7P2F "\n",
  238.         r6, ibeta, w);
  239.     printf("    OCCURRED FOR X =" F17P6E "\n", x1);
  240.     w = AMAX1(ait + w, ZERO);
  241.     printf(
  242.         " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS" F7P2F "\n\n\n",
  243.         ibeta, w);
  244.     w = -999.0e0;
  245.     if (r7 != ZERO)
  246.         w = ALOG(ABS(r7)) / albeta;
  247.     printf(" THE ROOT MEAN SQUARE RELATIVE ERROR WAS" F15P4E " = %4d **" F7P2F "\n",
  248.         r7, ibeta, w);
  249.     w = AMAX1(ait + w, ZERO);
  250.     printf(
  251.         " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS" F7P2F "\n\n\n",
  252.         ibeta, w);
  253.     a = 18.84955592e0L;
  254.     if (j == 2)
  255.         a = b + c;
  256.     b = a + c;
  257.     }
  258.  
  259.     /* special tests */
  260.  
  261.     printf("\fSPECIAL TESTS\n\n");
  262.     c = ONE / ipow(beta, (it / 2));
  263.     z = (sin(a + c) - sin(a - c)) / (c + c);
  264.     printf(" IF " F13P6E " IS NOT ALMOST 1.0E0,    SIN HAS THE WRONG PERIOD.\n\n",
  265.     z);
  266.  
  267.     printf(" THE IDENTITY   SIN(-X) = -SIN(X)   WILL BE TESTED.\n");
  268.     printf("     X        F(X) + F(-X)\n");
  269.  
  270.     for (i = 1; i <= 5; ++i)
  271.     {
  272.     x = ran(i1) * a;
  273.     z = sin(x) + sin(-x);
  274.     printf(F15P7E F15P7E "\n\n", x, z);
  275.     }
  276.  
  277.     printf(" THE IDENTITY SIN(X) = X , X SMALL, WILL BE TESTED.\n");
  278.     printf("    X         X - F(X)\n\n");
  279.     betap = ipow(beta, it);
  280.     x = ran(i1) / betap;
  281.  
  282.     for (i = 1; i <= 5; ++i)
  283.     {
  284.     z = x - sin(x);
  285.     printf(F15P7E F15P7E "\n\n", x, z);
  286.     x = x / beta;
  287.     }
  288.  
  289.     printf(" THE IDENTITY   COS(-X) = COS(X)   WILL BE TESTED.\n");
  290.     printf("        X         F(X) - F(-X)\n");
  291.  
  292.     for (i = 1; i <= 5; ++i)
  293.     {
  294.     x = ran(i1) * a;
  295.     z = cos(x) - cos(-x);
  296.     printf(F15P7E F15P7E "\n\n", x, z);
  297.     }
  298.  
  299.     printf(" TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT.\n\n");
  300.     expon = (float) minexp *0.75e0L;
  301.     x = pow(beta, expon);
  302.     y = sin(x);
  303.     printf("\n       SIN(" F13P6E ") = " F13P6E "\n", x, y);
  304.     printf(" THE FOLLOWING THREE LINES ILLUSTRATE THE LOSS IN SIGNIFICANCE\n");
  305.     printf(" FOR LARGE ARGUMENTS.  THE ARGUMENTS ARE CONSECUTIVE.\n");
  306.     z = sqrt(betap);
  307.     x = z * (ONE - epsneg);
  308.     y = sin(x);
  309.     printf("\n       SIN(" F13P6E ") =" F13P6E "\n", x, y);
  310.     y = sin(z);
  311.     printf("\n       SIN(" F13P6E ") =" F13P6E "\n", z, y);
  312.     x = z * (ONE + eps);
  313.     y = sin(x);
  314.     printf("\n       SIN(" F13P6E ") =" F13P6E "\n", x, y);
  315.  
  316.     /* test of error returns */
  317.  
  318.     printf("\fTEST OF ERROR RETURNS\n\n\n");
  319.  
  320.     x = betap;
  321.     printf(" SIN WILL BE CALLED WITH THE ARGUMENT" F15P4E "\n",x);
  322.     printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
  323.     fflush(stdout);
  324.     errno = 0;
  325.     y = sin(x);
  326.     if (errno)
  327.     perror("sin()");
  328.     printf(" SIN RETURNED THE VALUE" F15P4E "\n\n\n\n", y);
  329.  
  330.     printf(" THIS CONCLUDES THE TESTS\n");
  331. }
  332.